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ABSTRACT 

A new, high-order, conservative, and efficient method for conservation laws on unstructured 
grids is developed. It combines the best features of structured and unstructured grid methods to 
attain computational efficiency and geometric flexibility; it utilizes the concept of discontinuous 
and high-order local representations to achieve conservation and high accuracy; and it is based 
on the finiite-difference formulation for simplicity. Universal reconstructions are obtained by 
distributing unknowns in a geometrically similar manner for all unstructured cells. Placements 
of the unknown and flux points with various order of accuracy are given for the line, triangular 
and tetrahedral elements. The data structure of the new method permits an optimum use of cache 
memory, resulting in further computational efficiency on modern computers. A new pointer 
system is developed that reduces memory requirements and simplifies programming for any 
order of accuracy. Numerical solutions are presented and compared with the exact solutions for 
wave propagation problems in both two and three dimensions to demonstrate the capability of 
the method. Excellent agreement has been found. The method is simpler and more efficient than 
previous discontinuous Galerkin and spectral volume methods for unstructured grids. 


1. INTRODUCTION 

A current problem of great interest is to obtain very high-order accurate and efficient 
numerical solutions of systems of conservation laws for complex shapes. There are two types of 
computational grids, namely structured and unstructured grids. For most problems with complex 
shapes, unstructured grids consisting of triangles in 2D and tetrahedra in 3D are preferred 
because of the relatively short grid generation time. The most commonly used unstructured 
method is the finite-volume (FV) method, applied to the integral form of the conservation law 
with cell averages of the conservative variables as the unknowns. A reconstruction of any desired 
order of accuracy for each cell is obtained in terms of unknowns at neighboring cells. The flux 
integral for each face is evaluated using the reconstructed solution for the two cells and an 
approximate Riemann solver. For non-linear flux functions, a quadrature approximation is 
employed. Thus, conservation is satisfied locally for each cell. However, due to the unstructured 
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nature of the grid, it is difficult to obtain a non-singular stencil. This necessitates a least-squares 
inversion in general. For very high order of accuracy, the number of cells, and thus the number 
of operations to carry out the numerical procedure, can become very large in three dimensions. 
This would hamper the efficiency of the method. Furthermore, since each unknown employs a 
different stencil, one must repeat the least-squares inversion for every cell at each time step, or to 
store the inversion coefficients. In a high-order, three-dimensional computation, the former 
would involve impractically large CPU times, while for the latter the memory requirement 
becomes prohibitive. In addition, the data from neighboring cells required for the computation 
can be far apart in memory. This further hampers the efficiency of the method due to data 
gathering and scattering. As a result of these deficiencies, the FV method is limited to second 
order accuracy in most applications. 

Finite element (FE) methods have long been used for unstructured grids because of their 
geometric flexibility. A major difference between the FE and the FV is that the former seeks 
reconstruction data from within the cell (or element), while the latter employs the data from 
outside the cell. In the FE formulation, unknowns are nodal values, which are placed at 
geometrically similar points in each element. Due to this geometric similarity, the local 
reconstruction becomes universal for all cells in terms of the same set of cardinal basis or shape 
functions. By assigning a single value to the nodes on cell boundaries, a global piecewise 
continuous reconstruction can be obtained. The conservative equations are then satisfied in a 
weak form by multiplying them with the requisite number of test functions, integrating over the 
global domain, and using integration by parts. This results in a set of coupled equations for all 
the unknowns. Their solution involves a very large, sparse matrix, whose elements depend on the 
cell geometries. For nonlinear equations, quadrature approximations are necessary to evaluate 
the matrix elements. While the integral conservation law is satisfied for the global domain, it is 
not satisfied locally for each cell. Due to these deficiencies, the FE method is not as popular as 
the FV method for solving conservation laws. 

To overcome the deficiencies of the FE method for conservation laws, the discontinuous 
Galerkin (DG) method was introduced recently. Nodes on cell boundaries are allowed to have 
multiple values, so that the local reconstruction in each cell is in general discontinuous with that 
of its neighbors. The integration is carried out over each element, using the local shape functions 
as the test functions. As in the FV method, a Riemann solver is employed at cell boundaries. The 
integral conservation law is now satisfied for each cell. While we must still solve sets of coupled 
equations, each set involves now only the unknowns in one element. In addition, since the 
integrals involve quadratic terms, the quadrature formulas must have twice the degree of 
precision as the precision of the reconstruction. 

The spectral volume (SV) method overcomes the computational inefficiencies of the FV 
method by adopting the FE concept of using only data within the cell for reconstruction. Each 
triangular or tetrahedral cell, here called a spectral volume (SV), is partitioned into structured 
subcells called control volumes (CV). These are polygons in 2D, and polyhedra in 3D. The latter 
can have non-planar faces, which must be subdivided into planar facets in order to perform flux 
integrations. The unknowns are now cell averages over the CV's.If the SV's are partitioned in a 
geometrically similar manner, a single, universal reconstruction results. Thus only a few 
coefficients need to be stored in advance to evaluate all flux integrals. For high orders of 
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accuracy in 3D, the partitioning requires the introduction of a large number of parameters, whose 
optimization to achieve convergence becomes increasingly more difficult. The growth in the 
number of interior facets and the increase in the number of quadrature points for each facet, 
raises the computational cost greatly. In addition, the partitioning of curved cells, which are 
necessary at the boundaries representing complex shapes, is not straightforward. 

The finite difference (FD) method is based on the direct discretiazation of the differential 
form of the equations. While it is thus simple, it is not conservative, and for unstructured grids, a 
different reconstruction stencil is required at each point. In this respect it has the same 
disadvantages as the FV method. The advantages of the SV method and the simplicity of the FD 
method are combined in the new method to be presented below. 


2. THE SPECTRAL DIFFERENCE METHOD 

2. 1 General Description 

The partitioning of grid cells into subcells in the SV method is dictated by the need to 
satisfy the integral conservation law for each cell in order to capture discontinuities. It is then 
natural to define cell averages of conservative variables as discrete unknowns. But the flux 
integral for each grid cell can only be evaluated approximately. The accuracy of the flux at any 
point depends on the accuracy of the conservative variables, which is limited by the degree of 
precision of its reconstruction from the discrete unknowns. If the flux is non-linear, the surface 
integral is evaluated as a weighted sum of values at quadrature points, which is also accurate up 
to a certain degree of precision. Thus the unknowns can only be updated with an accuracy of a 
certain degree of precision. It is therefore sufficient to define the conservative unknowns at 
quadrature points that will approximate the volume integral over the cell to the desired order of 
accuracy, and that support a polynomial reconstruction of the desired degree of precision. 
Similarly, the flux will be defined at a different set of points, some or all of which lie on the cell 
boundary. The unknowns are updated using the differential form of the conservation law. 
Therefore the flux points must support a polynomial reconstruction of one higher degree of 
precision than the one for the conservative unknowns. Their location must also be chosen such 
that the integral conservation law for each grid cell is satisfied to the desired order of accuracy. If 
the nodes are distributed in a geometrically similar manner for all cells, the discretizations 
become universal, and can be expressed as the same weighted sums of the products of the local 
metrics and fluxes. These metrics are constants for the line, triangle, and tetrahedron elements, 
and can be computed analytically for curved elements. 


2.2 Details of the spectral difference method 


The most general form of a conservation law can be written as 

— + V*F = 0, 
dt 
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( 1 ) 


where the conservative variable u can be a scalar or a vector, and the generalized flux F can be a 
vector or tensor. The term V * F represents the divergence or curl of F, depending on the 
physical definition of u . Integrating (1) over cell i, we obtain 

— f udV+Jf dS*F = 0, (2) 

where V t is the volume of the cell i, and S u is the area of face l for cell i. (In 2D, each face is 
actually a line.) Here M is the number of faces, which is one more than the dimension D of the 
domain. 


In each cell, the discrete unknowns are the values of u at quadrature points for the 
volume integral over the cell. We denote these points, some of which may lie on the cell faces, as 
Tj ., and define . as u( r ; .) . The expansion of u in the cell can be written in the cardinal form 

u l (r) = ^L Ji (r)u JJ , (3) 

j - 1 

where L . .( r) are the cardinal basis functions and N n is the number of basis functions required to 

support a desired degree of precision n of the reconstruction. We will use polynomials as an 
independent basis. The locations of r y i then uniquely define the L ; ,.(r). 


In order to evaluate the surface integrals in (2) efficiently, we discretize F at points r ki , 
most or all of which are located at quadrature points for those integrals. If we define F ki as 

F(u(r u )), the expansion of F in the cell can also be written in the cardinal form 

w„„ 

k-1 

where M ki (r ) are now the set of cardinal basis functions defined by r k i . We can satisfy (1) at 
points i v . by evaluating 


• 'ii»i 

k - 1 

In order to evaluate F k i , u k i is required, which can be obtained directly from (3) 

u kj = ^ L jj( r u) u jj- 


(5) 

( 6 ) 
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To reduce the cost of this interpolation, some of points r k i may be chosen to coincide with r ;>j . If 
the points r . , and r k ; are distributed in a geometrically similar manner for all cells and within 

each cell i the gradient of a function is expressed in terms of its area vectors S-, the coefficients 
in (5) and *(6) become universal, independent of cell i. Therefore, one can express these 
coefficients as 


and 


Lj,iir k .) = l kj , 

V i 1-1 


(7) 

( 8 ) 
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There are only a few of these coefficients, which can be calculated exactly and stored in advance. 
Equations (6) and (5) can thus be rewritten as 


U k,i ~ 2 hj U j,i ’ 


(9) 


j - 1 


and 


V * f (r J.<> ~ ^ 2 S .' * 2 (10) 

V i l-\ Jc-l 

For those points r k , located on the cell faces, since u may be discontinuous, we must replace the 
flux F by a Riemann flux. Let r tv . be the corresponding F point in the neighboring cell i . 
Define u L as u,.( r ki ), and u R as u,.,(r tv .), and let F L =F(u L ) and F R = F(u R ). We further define 

d(n * F ) 


n as the unit vector outward normal to the face and A ; 


Riemann flux is then given by 


du 


as the Jacobian matrix. The 


Fjliem ~ ^ [^R + ^ L~ ^ ^ ^ ( U R U L ) H ] » 

where A(u L ,u R ,n) is defined by the particular Riemann solver. 

In order to check the integral conservation law for the cell, we evaluate 

N n 

f v V*FdV = V i 2w J V*F(r J , i ) 
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using Eq. (10). Here w . are the volume quadrature weights at the points iv,-. We can also 
evaluate the surface integral of F as 

M M 

< s *' ? -2 s !*2 b v f *-,- < I3 > 

l-\ k 

(9) 

Here w k are surface quadrature weights for face l, and for each face the summation is over those 

D 

points r tii located on that face. Note that Sf = - 2*: . Conservation requires that (24) and (25) 

N! 

yield the same value, within the accuracy of the method. The total contribution from any interior 
point F k i to the volume integral in (12) must vanish. 

The SD formulation for the line element is a generalization of the multidomain spectral 
method. Here the distributions of u and F points are based on those quadrature points that lead 
to the satisfaction of the integral form of the conservation law. Since this does not define a 
unique set of locations, they can be optimized by minimizing some object functions. The tensor 
products of the line formulation can be used for quadrilateral and hexahedral cells or elements. 



From Eqs. (9) and (10), we also see that the SD formulation is very similar to that of the 
FD method for structured grids. The SD method thus retains the simplicity and computational 
efficiency of the structured FD method. However, the metric terms in the latter are evaluated by 
numerically differencing the grid point coordinates. Since numerical grid generators are mostly 
only second-order accurate, the overall accuracy of the solution can be severely degraded if the 
grid is not sufficiently smooth. In contrast, the metric terms in the SD method are computed 
exactly from the geometry of the grid, no matter how it was generated. It thus retains its formal 
accuracy, even for very unsmooth unstructured grids. Furthermore, in contrast to the FD method, 
the integral conservation law is satisfied to the desired accuracy. 

2.3 Data Structure 

There are several aspects of the data structure which lead to a very efficient parallelizable 
code. The code utilizes a new pointer system that reduces memory requirements and simplifies 
the programming for any order of accuracy. The SD method permits an optimum use of Cache 
memory, since all the data required to perform the calculations for a given cell are found 
contiguous in memory, resulting in great reductions in memory access. Both of these subjects 
will be discussed in more details in the final paper. 


3. LOCATIONS OF THE UNKNOWN AND FLUX POINTS 

The critical part of the SD method is the location of the u points r jt and F points v ki . 

If those points are distributed in a geometrically similar manner for all cells, the formulas for the 
flux derivatives become universal, and can be expressed as the same weighted sums of the 
products of the local metrics and fluxes. In the final paper, we shall indicate how these metrics 
are evaluated for a general curvilinear tetrahedron. (A curvilinear triangle is just a special case). 
All the relations for cells with planar faces are then valid in terms of the transformed u and F 
values at quadrature points. 

The locations of the u and F points are determined by symmetry groups associated with 
the cell centroids, vertices, edges, and faces. All but the first contain arbitrary parameters that 
can be varied to obtain optimum solutions. The number of points required for a reconstruction 
with a specified degree of precision is greater than the minimum number of Gaussian quadrature 
points for that precision. One can obtain greater efficiency by locating some u points to coincide 
with F points. For F points on the vertices in 2D and edges and vertices in 3D, more than one 
Riemann solver is necessary. For those formulations with expensive Riemann solvers, these 
points should be minimized. Another criterion for the placement of u and F points is that the 
reconstruction matrix is nonsingular. The final criterion is that integral conservation is satisfied 
within the desired degree of precision. In order to achieve a nonsingular matrix, one could allow 
an overspecifying of the number of F points. The reconstruction would then require a least- 
squares inversion. Since this occurs only in a preprocessing phase, it would not hamper the 
efficiency of the method. Even in this case, we can show that the number of F points is far less 
than the number of flux quadrature points in the SV method with the same accuracy. 
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We first show some representative placements of u points (circles) and F points 
(squares) that satisfy integral conservation, for various orders of accuracy (which are one more 
than the degree of precision). The placements for the line element are shown in Fig. la for 
second order and Fig. lb for fifth order accuracy. The u points are located at the Gaussian 
quadrature points. Note that the same order of accuracy is achieved by locating the u points at 
Gauss-Lobatto quadrature points. While we now have one more u unknown per element, we 
eliminate all interpolations by locating all u points at the F points. 

In Fig. 2 we show the placement of the points for the triangle. Fig. 2a and 2b show two 
different placements of the F points for first order accuracy. Three different second order 
accurate placements are shown. In Fig. 2c, the u points are at the Gaussian quadrature points, 
while the F points are at the Gauss-Lobatto points on the edges. In Fig. 2d the u points are 
moved to the other Gauss points on the midedges, where they coincide with F points. Moving 
the F points to the Gaussian points on the edges in order to minimize Riemann solver calls 
results in a singular matrix. The singularity is removed by adding an additional F point at the 
centroid. This case is shown in Fig. 2e. Two third order accurate placements are shown. In Fig. 
2f, the u points are at the Gaussian points, while the F points are at the Gauss-Lobatto points on 
the edges. In Fig. 2g the F points on the edges have been moved to the Gaussian quadrature 
points, to reduce the number of Riemann solver calculations. 

In Fig. 3 we show several placements of the points for the tetrahedron. Fig. 3a and 3b 
show two different placements of the F points for first order accuracy. A second order accurate 
placement of points is shown in Fig. 3c. The u points are located at the Gaussian quadrature 
points. Six of the F points are on the edge midpoints, while the remaining four vertex based 
points are at some arbitrary interior points. 

In the final paper we will show the results of optimization studies and define the 
optimized placements of the u and F points for line, triangles and tetrahedra for various orders 
of accuracy. 


4. NUMERICAL RESULTS 

We first show some one-dimensional convergence studies for the linear convection 
equation. Plotted are L, (solid lines) and L x (dash lines) error norms as functions of grid size for 
second to fifth order accurate methods. Fig. 4a shows the plots for the u points at the Gaussian 
quadrature points, while Fig. 4b shows the plots for the u points at the Gauss-Lobatto quadrature 
points. Both exhibit the expected orders of accuracy. In Fig. 5 we show analogous plots for a 
plane wave propagating at 45 degree through a square domain for first to third order accurate 
methods, using the point placements in Figs. 2a, 2c, and 2f. The third example is the scattering of 
a TM wave incident on a perfectly conducting circular cylinder. Fig. 6 shows the unstructured 
grid consisting of 2024 triangular cells. The wave propagates from the left to the right with the 
wave number equal to 5, based on the radius of the cylinder. This gives approximately 6 cells per 
wavelength. Contour plots for E z with the exact solution (solid lines) are shown in Figs. 7a, 7b, 
and 7c for first, second, and third order accurate methods, respectively. It is seen that the first 
order solution is very dissipative with this grid resolution, while the second and third order 
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solutions show an excellent agreement with the exact solution. In the final paper we will show 
additional high-order 2D as well as 3D results, including curved boundaries. 

5. CONCLUDING REMARKS 

We presented a very efficient finite-difference method for unstructured grids that satisfies 
conservation for high orders of accuracy. The structured distribution of discrete variables in each 
unstructured cell maintains computational efficiency and geometric flexibility. The flux 
derivatives needed to update the conservative unknowns are expressed as universal weighted 
sums of the unknowns, leading to great computational efficiency. An important aspect of the 
method is that the number of Riemann solvers per unknown decreases as the order of accuracy 
increases, reducing the cost for higher order. Curved boundaries can be implemented in a simple 
manner. The data structure permits an optimum use of cache memory, resulting in great 
reduction in data communication. A new pointer system is used that reduces memory 
requirements and simplifies programming for any order of accuracy. The method is simpler and 
more efficient than previous DG and SV methods for unstructured grids. Further improvements 
include implicit time integration, grid and order adaptation, multidimensional limiters and filters, 
and moving grids. 
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Figure 1. Placement of unknown and flux points for a line element 
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Figure 3. Placement 



Figure 5. Error norms 








Figure 7. Contour plots of E z for a plane wave incident on a perfectly conducting cylinder 
(numerical solutions: color contours, exact solution: solid lines) 


(a) first order 


Figure 6. Grid for a region exterior to a circular cylinder. 

(2) second order (c) third order 





